# analysis replication for Dawson, Schwenk and Xezonakis (2024). The distribution of executive power and corruption: A meta-analytical review
rm(list=ls())

# libraries                 # version used in this analysis
library(psych)              # 2.3.9
library(tidyverse)          # 2.0.0
library(ggpubr)             # 0.6.0
library(metafor)            # 4.6-0
library(broom)              # 1.0.5
library(dotwhisker)         # 0.8.1
library(cowplot)            # 1.1.3

# set working directory

# load data
stdat <- read_csv("studies_all.csv") # STUDY-LEVEL
edat.all <- read_csv("estimates_all.csv") # ESTIAMTE-LEVEL - ALL

# Subsets for each IV
edat.pres <- edat.all %>% 
  filter(ivtype=="presidentialism")

edat.decpol <- edat.all %>% 
  filter(ivtype=="decent_pol")

edat.decfis <- edat.all %>% 
  filter(ivtype=="decent_fisc")

################################################################################

############################## FIGURE 1 ########################################

length(unique(edat.all$idstud))
length(unique(edat.all$idmod))
length(edat.all$idest)

## PANEL A
describe(stdat$yrstud)
fig1_a <- stdat %>% 
  ggplot(aes(x = yrstud)) +
  geom_bar(color = "black", fill = "#939593") +
  labs(x = "", 
       y = "", 
       title = "A. Year of publication") +
  scale_y_continuous(breaks = c(0,2,4,6),
                     limits = c(0,6)) +
  scale_x_continuous(breaks = c(2000,2005,2010,2015,2020,2025)) +
  theme_bw() +
  theme(plot.title = element_text(face = "bold"),
        text = element_text(size = 15))

## PANEL B
table(edat.all$dv_measure)
fig1_b <- edat.all %>% 
  mutate(dv_measure = factor(dv_measure)) %>% 
  count(dv_measure) %>% 
  mutate(perc = n/nrow(edat.all)) %>% 
  ggplot(aes(x = reorder(dv_measure, -perc), y = perc)) +
  geom_bar(stat = "identity", color = "black", fill = "#939593") +
  labs(x = "", 
       y = "",
       title = "B. Corruption measure") +
  scale_y_continuous(labels = scales::percent, 
                     limits = c(0,.4),
                     breaks = c(0,.1,.2,.3,.4)) +
  theme_bw() +
  theme(plot.title = element_text(face = "bold"),
        text = element_text(size = 15))

## PANEL C
fig1_c <- edat.all %>% 
  group_by(autyr) %>% 
  summarise_at(vars(yrstart, yrend, ds.ts), 
               list(group_yr = mean, mode = function(x){
                 as.numeric(names(sort(-table(x)))[1])})) %>%
  mutate(autyr = factor(autyr)) %>% 
  mutate(ds.ts_mode = factor(ds.ts_mode)) %>% 
  mutate(autyr = fct_reorder(autyr, yrstart_group_yr)) %>% 
  ggplot() +
  geom_linerange(aes(x = autyr, 
                     ymin = yrstart_group_yr, 
                     ymax = yrend_group_yr,
                     color = ds.ts_mode),
                 size = 1) +
  geom_point(aes(x = autyr,
                 y = yrstart_group_yr, color = ds.ts_mode), 
             shape = 15, size = 1) +
  geom_point(aes(x = autyr,
                 y = yrend_group_yr, color = ds.ts_mode), 
             shape = 15, size = 1) +
  labs(x = "Study", 
       y = "", 
       color = "", 
       title = "C. Data coverage across studies") +
  scale_y_continuous(breaks = c(1960,1980,2000,2020)) +
  scale_color_manual(labels = c("Cross section", "Time series"), 
                     values = c("0" = "black", "1"="grey")) +
  theme_bw() +
  theme(axis.text.x = element_blank(),
        plot.title = element_text(face = "bold"),
        legend.position = "bottom",
        text = element_text(size = 15))

fig1 <- ggarrange(ggarrange(fig1_a, fig1_b, ncol = 2),
                  fig1_c, nrow = 2)

################################################################################

############ META-ANALYSIS - OVERALL (APPENDIX C, TABLE 5) #####################

# Presidentialism
m.pres_base <- rma.mv(yi = partcorr, 
                      V = var_partcorr,
                      slab = idstud,
                      data = edat.pres,
                      random = ~ 1 | idstud/idest, 
                      test = "t")
summary(m.pres_base)

# Political Decentralization
m.decpol_base <- rma.mv(yi = partcorr, 
                        V = var_partcorr,
                        slab = idstud,
                        data = edat.decpol,
                        random = ~ 1 | idstud/idest, 
                        test = "t")
summary(m.decpol_base)

# Fiscal decentralization
m.decfis_base <- rma.mv(yi = partcorr, 
                        V = var_partcorr,
                        slab = idstud,
                        data = edat.decfis,
                        random = ~ 1 | idstud/idest, 
                        test = "t")
summary(m.decfis_base)

############ META-ANALYSIS - FOCAL VARS (APPENDIX C, TABLE 6) ##################

# Presidentialism #
m.pres_nocontrols <- rma.mv(yi = partcorr, 
                            V = var_partcorr,
                            slab = idstud,
                            data = edat.pres %>% filter(control.no==1),
                            random = ~ 1 | idstud/idest, 
                            test = "t")
summary(m.pres_nocontrols)

# Political Decentralization #
m.decpol_nocontrols <- rma.mv(yi = partcorr, 
                              V = var_partcorr,
                              slab = idstud,
                              data = edat.decpol %>% filter(control.no==1),
                              random = ~ 1 | idstud/idest, 
                              test = "t")
summary(m.decpol_nocontrols)

# Fiscal decentralization #
m.decfis_nocontrols <- rma.mv(yi = partcorr, 
                              V = var_partcorr,
                              slab = idstud,
                              data = edat.decfis %>% filter(control.no==1),
                              random = ~ 1 | idstud/idest, 
                              test = "t")
summary(m.decfis_nocontrols)


################################################################################

############################## FIGURE 2 ########################################

fig2_a <- edat.pres %>% 
  arrange(partcorr) %>% 
  mutate(id = seq.int(nrow(edat.pres))) %>% 
  ggplot() +
  geom_hline(yintercept = 0) +
  geom_hline(yintercept = 0.08, color = "red") +
  annotate("text", x = 25, y = 0.6,
           label = "b = 0.08") +
  annotate("text", x = 17, y = 0.6,
           label = "k = 163") +
  annotate("text", x = 9, y = 0.6,
           label = "t = 2.56") +
  geom_hline(yintercept = 0.018, linetype='dotted', color = "red") +
  geom_hline(yintercept = 0.138, linetype='dotted', color = "red") +
  geom_linerange(aes(x = id,
                     ymin = ci_low,
                     ymax = ci_upp),
                 color = "black", alpha = 0.3) +
  geom_point(aes(x = id, y = partcorr)) +
  labs(x = "Estimate", y = "", title = "Presidentialism") +
  scale_y_continuous(limits=c(-1, 1)) +
  coord_flip() +
  theme_bw()  +
  theme(axis.text.y=element_blank(),
        plot.title = element_text(face = "bold"),
        text = element_text(size = 15))

fig2_b <- edat.decpol %>% 
  arrange(partcorr) %>% 
  mutate(id = seq.int(nrow(edat.decpol))) %>% 
  ggplot() +
  geom_hline(yintercept = 0) +
  geom_hline(yintercept = 0.063, color = "red") +
  annotate("text", x = 63, y = 0.6,
           label = "b = 0.063") +
  annotate("text", x = 42, y = 0.6,
           label = "k = 470") +
  annotate("text", x = 21, y = 0.6,
           label = "t = 2.49") +
  geom_hline(yintercept = 0.013, linetype='dotted', color = "red") +
  geom_hline(yintercept = 0.112, linetype='dotted', color = "red") +
  geom_linerange(aes(x = id,
                     ymin = ci_low,
                     ymax = ci_upp),
                 color = "black", alpha = 0.3) +
  geom_point(aes(x = id, y = partcorr)) +
  labs(x = "", y = "Partial correlation", title = "Political Decentralization") +
  scale_y_continuous(limits=c(-1, 1)) +
  coord_flip() +
  theme_bw()  +
  theme(axis.text.y=element_blank(),
        plot.title = element_text(face = "bold"),
        text = element_text(size = 15))


fig2_c <- edat.decfis %>% 
  arrange(partcorr) %>% 
  mutate(id = seq.int(nrow(edat.decfis))) %>% 
  ggplot() +
  geom_hline(yintercept = 0) +
  geom_hline(yintercept = -0.116, color = "red") +
  annotate("text", x = 50, y = 0.6,
           label = "b = -0.116") +
  annotate("text", x = 34, y = 0.6,
           label = "k = 379") +
  annotate("text", x = 18, y = 0.6,
           label = "t = -3.91") +
  geom_hline(yintercept = -0.175, linetype='dotted', color = "red") +
  geom_hline(yintercept = -0.058, linetype='dotted', color = "red") +
  geom_linerange(aes(x = id,
                     ymin = ci_low,
                     ymax = ci_upp),
                 color = "black", alpha = 0.3) +
  geom_point(aes(x = id, y = partcorr)) +
  labs(x = "", y = "", title = "Fiscal Decentralization") +
  scale_y_continuous(limits=c(-1, 1)) +
  coord_flip() +
  theme_bw()  +
  theme(axis.text.y=element_blank(),
        plot.title = element_text(face = "bold"),
        text = element_text(size = 15))

fig2 <- ggarrange(fig2_a, fig2_b, fig2_c,
                  ncol = 3, nrow = 1, align = "h")

################################################################################

############################## FIGURE 3 ########################################

################ NUMERICAL RESULTS FOR APPENDIX C, TABLES 7-9 ##################

## COLUMN 1: PRESIDENTIALISM ##
## ROW A: DATA STRUCTURE
m.pres_ds <- rma.mv(yi = partcorr, 
                    V = var_partcorr, 
                    mods = cbind(ds.cs99, ds.cs00, ds.ts),
                    slab = idstud,
                    data = edat.pres,
                    random = ~ 1 | idstud/idest, 
                    intercept = FALSE,
                    test = "t")
summary(m.pres_ds)
anova(m.pres_base, m.pres_ds, refit = TRUE)
fig3_1a <- m.pres_ds %>% 
  tidy() %>%
  relabel_predictors(c(ds.cs99 = "XS <2000",
                       ds.cs00 = "XS ≥2000",
                       ds.ts = "Time series")) %>% 
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7)) +
  annotate("rect", xmin = 0.018, xmax = 0.138, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "Presidentialism") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  theme(axis.text.x=element_blank(),
        legend.position = "none",
        plot.margin=grid::unit(c(0,1,0,0), "mm"),
        plot.title = element_text(face = "bold"),
        text = element_text(size = 15))

## ROW B: METHODOLOGY
m.pres_meth <- rma.mv(yi = partcorr, 
                      V = var_partcorr, 
                      mods = cbind(meth.ols, meth.wls, meth.2sls, meth.ord, meth.oth),
                      slab = idstud,
                      data = edat.pres,
                      random = ~ 1 | idstud/idest,
                      intercept = FALSE,
                      test = "t")
summary(m.pres_meth)
anova(m.pres_base, m.pres_meth, refit = TRUE)
fig3_1b <- m.pres_meth %>% 
  tidy() %>%
  add_row(term = "2SLS", type = NA, estimate = NA, std.error = NA, statistic = NA, p.value = NA) %>% 
  relabel_predictors(c(meth.ols = "OLS",
                       meth.wls = "WLS",
                       meth.ord = "Ordered",
                       meth.oth = "Other")) %>% 
  arrange(factor(term, levels = c("OLS", "WLS", "2SLS", "Ordered", "Other"))) %>% 
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7)) +
  annotate("rect", xmin = 0.018, xmax = 0.138, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  theme(axis.text.x=element_blank(),
        legend.position = "none",
        plot.margin=grid::unit(c(-5,1,0,0), "mm"),
        text = element_text(size = 15))

## ROW C: SAMPLE
m.pres_sam <- rma.mv(yi = partcorr, 
                     V = var_partcorr, 
                     mods = cbind(sam_global, sam_developed, sam_developing),
                     slab = idstud,
                     data = edat.pres,
                     random = ~ 1 | idstud/idest,
                     intercept = FALSE,
                     test = "t")
summary(m.pres_sam)
anova(m.pres_base, m.pres_sam, refit = TRUE)
fig3_1c <- m.pres_sam %>% 
  tidy() %>%
  relabel_predictors(c(sam_global = "Global",
                       sam_developed = "Developed",
                       sam_developing = "Developing")) %>% 
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7)) +
  annotate("rect", xmin = 0.018, xmax = 0.138, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  theme(legend.position = "none",
        plot.margin=grid::unit(c(-5,1,0,0), "mm"),
        text = element_text(size = 15))

## COLUMN 2: POLITICAL DECENTRALIZATION
## ROW A: DATA STRUCTURE
m.decpol_ds <- rma.mv(yi = partcorr, 
                      V = var_partcorr, 
                      mods = cbind(ds.cs99, ds.cs00, ds.ts),
                      slab = idstud,
                      data = edat.decpol,
                      random = ~ 1 | idstud/idest,
                      intercept = FALSE,
                      test = "t")
summary(m.decpol_ds)
anova(m.decpol_base, m.decpol_ds, refit = TRUE)
fig3_2a <- m.decpol_ds %>% 
  tidy() %>%
  relabel_predictors(c(ds.cs99 = "XS <2000",
                       ds.cs00 = "XS ≥2000",
                       ds.ts = "Time series")) %>% 
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7)) +
  annotate("rect", xmin = 0.013, xmax = 0.112, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "Pol. Decentralization") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        legend.position = "none",
        plot.margin=grid::unit(c(0,1,0,0), "mm"),
        plot.title = element_text(face = "bold"),
        text = element_text(size = 15))

## ROW B: METHODOLOGY
m.decpol_meth <- rma.mv(yi = partcorr, 
                        V = var_partcorr, 
                        mods = cbind(meth.ols, meth.wls, meth.2sls, meth.ord, meth.oth),
                        slab = idstud,
                        data = edat.decpol,
                        random = ~ 1 | idstud/idest, 
                        intercept = FALSE,
                        test = "t")
summary(m.decpol_meth)
anova(m.decpol_base, m.decpol_meth, refit = TRUE)
fig3_2b <- m.decpol_meth %>% 
  tidy() %>%
  relabel_predictors(c(meth.ols = "OLS",
                       meth.wls = "WLS",
                       meth.2sls = "2SLS",
                       meth.ord = "Ordered",
                       meth.oth = "Other")) %>% 
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7)) +
  annotate("rect", xmin = 0.013, xmax = 0.112, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        legend.position = "none",
        plot.margin=grid::unit(c(-5,1,0,0), "mm"),
        text = element_text(size = 15))

## ROW C: SAMPLE
m.decpol_sam <- rma.mv(yi = partcorr, 
                       V = var_partcorr, 
                       mods = cbind(sam_global, sam_developed, sam_developing),
                       slab = idstud,
                       data = edat.decpol,
                       random = ~ 1 | idstud/idest,
                       intercept = FALSE,
                       test = "t")
summary(m.decpol_sam)
anova(m.decpol_base, m.decpol_sam, refit = TRUE)
fig3_2c <- m.decpol_sam %>% 
  tidy() %>%
  relabel_predictors(c(sam_global = "Global",
                       sam_developed = "Developed",
                       sam_developing = "Developing")) %>% 
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7)) +
  annotate("rect", xmin = 0.013, xmax = 0.112, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  theme(axis.text.y=element_blank(),
        legend.position = "none",
        plot.margin=grid::unit(c(-5,1,0,0), "mm"),
        text = element_text(size = 15))

## COLUMN 3: FISCAL DECENTRALIZATION
## ROW A: DATA STRUCTURE
m.decfis_ds <- rma.mv(yi = partcorr, 
                      V = var_partcorr, 
                      mods = cbind(ds.cs99, ds.cs00, ds.ts),
                      slab = idstud,
                      data = edat.decfis,
                      random = ~ 1 | idstud/idest,
                      intercept = FALSE,
                      test = "t")
summary(m.decfis_ds)
anova(m.decfis_base, m.decfis_ds, refit = TRUE)
fig3_3a <- m.decfis_ds %>% 
  tidy() %>%
  relabel_predictors(c(ds.cs99 = "XS <2000",
                       ds.cs00 = "XS ≥2000",
                       ds.ts = "Time series")) %>% 
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7)) +
  annotate("rect", xmin = -0.175, xmax = -0.038, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "Fisc. Decentralization") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        legend.position = "none",
        plot.margin=grid::unit(c(0,1,0,0), "mm"),
        plot.title = element_text(face = "bold"),
        text = element_text(size = 15)) +
  ylab("Data Structure") +
  scale_y_discrete(position = "right")

## ROW B: METHODOLOGY
m.decfis_meth <- rma.mv(yi = partcorr, 
                        V = var_partcorr, 
                        mods = cbind(meth.ols, meth.wls, meth.2sls, meth.ord, meth.oth),
                        slab = idstud,
                        data = edat.decfis,
                        random = ~ 1 | idstud/idest, 
                        intercept = FALSE,
                        test = "t")
summary(m.decfis_meth)
anova(m.decfis_base, m.decfis_meth, refit = TRUE)
fig3_3b <- m.decfis_meth %>% 
  tidy() %>%
  relabel_predictors(c(meth.ols = "OLS",
                       meth.wls = "WLS",
                       meth.2sls = "2SLS",
                       meth.ord = "Ordered",
                       meth.oth = "Other")) %>% 
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7)) +
  annotate("rect", xmin = -0.175, xmax = -0.038, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        legend.position = "none",
        plot.margin=grid::unit(c(-5,1,0,0), "mm"),
        text = element_text(size = 15)) +
  ylab("Method") +
  scale_y_discrete(position = "right")

## ROW C: SAMPLE
m.decfis_sam <- rma.mv(yi = partcorr, 
                       V = var_partcorr, 
                       mods = cbind(sam_global, sam_developed, sam_developing),
                       slab = idstud,
                       data = edat.decfis,
                       random = ~ 1 | idstud/idest, 
                       intercept = FALSE,
                       test = "t")
summary(m.decfis_sam)
anova(m.decfis_base, m.decfis_sam, refit = TRUE)
fig3_3c <- m.decfis_sam %>% 
  tidy() %>%
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7)) +
  annotate("rect", xmin = -0.175, xmax = -0.038, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  theme(axis.text.y=element_blank(),
        legend.position = "none",
        plot.margin=grid::unit(c(-5,1,0,0), "mm"),
        text = element_text(size = 15)) +
  ylab("Sample") +
  scale_y_discrete(position = "right")

fig3 <- ggarrange(ggarrange(fig3_1a, fig3_1b, fig3_1c, nrow = 3, align = "v"),
                  NULL,
                  ggarrange(fig3_2a, fig3_2b, fig3_2c, nrow = 3),
                  NULL,
                  ggarrange(fig3_3a, fig3_3b, fig3_3c, nrow = 3), 
                  ncol = 5, widths = c(1, .05, .8, .1, .8))
fig3 <- ggdraw(fig3) +
  theme(plot.background = element_rect(fill="white", color = NA))

################################################################################

############################## FIGURE 4 ########################################

################ NUMERICAL RESULTS FOR APPENDIX C, TABLES 7-9 ##################

## COLUMN 1: PRESIDENTIALISM
## ROW A: IV MEASURE
m.pres_iv <- rma.mv(yi = partcorr, 
                    V = var_partcorr, 
                    mods = cbind(iv_pres_dummy, iv_pres_cat),
                    slab = idstud,
                    data = edat.pres,
                    random = ~ 1 | idstud/idest, 
                    intercept = FALSE,
                    test = "t")
summary(m.pres_iv)
anova(m.pres_base, m.pres_iv, refit = TRUE)
fig4_1a <- m.pres_iv %>% 
  tidy() %>%
  relabel_predictors(c(iv_pres_dummy = "Dummy",
                       iv_pres_cat = "Categorical")) %>% 
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7)) +
  annotate("rect", xmin = 0.018, xmax = 0.138, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "Presidentialism") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  theme(axis.text.x=element_blank(),
        legend.position = "none",
        plot.margin=grid::unit(c(0,1,0,0), "mm"),
        plot.title = element_text(face = "bold", hjust = 1),
        text = element_text(size = 15))

## ROW B: DV MEASURE
m.pres_dv <- rma.mv(yi = partcorr, 
                    V = var_partcorr, 
                    mods = cbind(dv.icrg, dv.wbank, dv.ti, dv.oth),
                    slab = idstud,
                    data = edat.pres,
                    random = ~ 1 | idstud/idest, 
                    intercept = FALSE,
                    test = "t")  
summary(m.pres_dv)
anova(m.pres_base, m.pres_dv, refit = TRUE)
fig4_1b <- m.pres_dv %>% 
  tidy() %>%
  relabel_predictors(c(dv.icrg = "ICRG",
                       dv.wbank = "World Bank",
                       dv.ti = "CPI",
                       dv.oth = "Other")) %>% 
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7)) +
  annotate("rect", xmin = 0.018, xmax = 0.138, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  theme(legend.position = "none",
        plot.margin=grid::unit(c(-5,1,0,0), "mm"),
        text = element_text(size = 15))

## COLUMN 2: POLITICAL DECENTRALIZATION
## ROW A: IV MEASURE
m.decpol_base_comp <- rma.mv(yi = partcorr, 
                             V = var_partcorr, 
                             slab = idstud,
                             data = edat.decpol %>% filter(!is.na(iv_decpol_dummy)),
                             random = ~ 1 | idstud/idest,
                             method = "ML")
m.decpol_iv <- rma.mv(yi = partcorr, 
                      V = var_partcorr, 
                      mods = cbind(iv_decpol_index, iv_decpol_dummy, iv_decpol_unit, iv_decpol_oth),
                      slab = idstud,
                      data = edat.decpol %>% filter(!is.na(iv_decpol_dummy)),
                      random = ~ 1 | idstud/idest, 
                      intercept = FALSE,
                      method = "ML")
summary(m.decpol_iv)
anova(m.decpol_base_comp, m.decpol_iv, refit = TRUE)
fig4_2a <- m.decpol_iv %>% 
  tidy() %>%
  relabel_predictors(c(iv_decpol_dummy = "Dummy",
                       iv_decpol_index = "Index",
                       iv_decpol_unit = "No. Units",
                       iv_decpol_oth = "Other")) %>% 
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7)) +
  annotate("rect", xmin = 0.013, xmax = 0.112, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "Pol. Decentralization") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  theme(axis.text.x=element_blank(),
        legend.position = "none",
        plot.margin=grid::unit(c(0,1,0,0), "mm"),
        plot.title = element_text(face = "bold", hjust = 1),
        text = element_text(size = 15))

## ROW B: DV MEASURE
m.decpol_dv <- rma.mv(yi = partcorr, 
                      V = var_partcorr, 
                      mods = cbind(dv.icrg, dv.wbank, dv.ti, dv.oth),
                      slab = idstud,
                      data = edat.decpol,
                      random = ~ 1 | idstud/idest, 
                      intercept = FALSE,
                      test = "t")
summary(m.decpol_dv)
anova(m.decpol_base, m.decpol_dv, refit = TRUE)
fig4_2b <- m.decpol_dv %>% 
  tidy() %>%
  relabel_predictors(c(dv.icrg = "ICRG",
                       dv.wbank = "World Bank",
                       dv.ti = "CPI",
                       dv.oth = "Other")) %>% 
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7)) +
  annotate("rect", xmin = 0.013, xmax = 0.112, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  theme(legend.position = "none",
        plot.margin=grid::unit(c(-5,1,0,0), "mm"),
        text = element_text(size = 15))

## COLUMN 3: FISCAL DECENTRALIZATION
## ROW A: IV MEASURE 
m.decfis_base_comp <- rma.mv(yi = partcorr, 
                             V = var_partcorr, 
                             slab = idstud,
                             data = edat.decfis %>% filter(!is.na(iv_decfis_oth)),
                             random = ~ 1 | idstud/idest,
                             method = "ML")
m.decfis_iv <- rma.mv(yi = partcorr, 
                      V = var_partcorr, 
                      mods = cbind(iv_decfis_rev, iv_decfis_exp, iv_decfis_index, iv_decfis_oth),
                      slab = idstud,
                      data = edat.decfis %>% filter(!is.na(iv_decfis_oth)),
                      random = ~ 1 | idstud/idest, 
                      intercept = FALSE,
                      method = "ML")
summary(m.decfis_iv)
anova(m.decfis_base_comp, m.decfis_iv, refit = TRUE)
fig4_3a <- m.decfis_iv %>% 
  tidy() %>%
  relabel_predictors(c(iv_decfis_rev = "Revenue",
                       iv_decfis_exp = "Expenditure",
                       iv_decfis_index = "Index",
                       iv_decfis_oth = "Other")) %>% 
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7)) +
  annotate("rect", xmin = -0.175, xmax = -0.038, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "Fisc. Decentralization") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  labs(tag = "IV Measure") +
  theme(axis.text.x=element_blank(),
        legend.position = "none",
        plot.margin=grid::unit(c(0,15,0,0), "mm"),
        plot.title = element_text(face = "bold", hjust = 1),
        text = element_text(size = 15),
        plot.tag=element_text(size = 13, angle=-90),
        plot.tag.position=c(1.05, 0.5))

## ROW B: DV MEASURE
m.decfis_dv <- rma.mv(yi = partcorr, 
                      V = var_partcorr, 
                      mods = cbind(dv.icrg, dv.wbank, dv.ti, dv.oth),
                      slab = idstud,
                      data = edat.decfis,
                      random = ~ 1 | idstud/idest, 
                      intercept = FALSE,
                      test = "t")
summary(m.decfis_dv)
anova(m.decfis_base, m.decfis_dv, refit = TRUE)
fig4_3b <- m.decfis_dv %>% 
  tidy() %>%
  relabel_predictors(c(dv.icrg = "ICRG",
                       dv.wbank = "World Bank",
                       dv.ti = "CPI",
                       dv.oth = "Other")) %>% 
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7)) +
  annotate("rect", xmin = -0.175, xmax = -0.038, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  labs(tag = "DV Measure") +
  theme(legend.position = "none",
        plot.margin=grid::unit(c(-5,15,0,0), "mm"),
        plot.title = element_text(face = "bold", hjust = 1),
        text = element_text(size = 15),
        plot.tag=element_text(size = 13, angle=-90),
        plot.tag.position=c(1.05, 0.5))

#combine figs
fig4 <- ggarrange(fig4_1a, fig4_2a, fig4_3a,
                  fig4_1b, fig4_2b, fig4_3b,
                  ncol = 3,
                  nrow = 2,
                  align = "v")

################################################################################

############################## FIGURE 5 ########################################

################ NUMERICAL RESULTS FOR APPENDIX C, TABLES 7-9 ##################

## COLUMN 1: PRESIDENTIALISM
m.pres_con <- rma.mv(yi = partcorr, 
                     V = var_partcorr, 
                     mods = cbind(con.bstruc, con.prop, con.distm, con.decpol, con.decfis),
                     slab = idstud,
                     data = edat.pres,
                     random = ~ 1 | idstud/idest, 
                     intercept = FALSE,
                     test = "t",
                     control=list(rel.tol=1e-8))
summary(m.pres_con)
anova(m.pres_base, m.pres_con, refit = TRUE)
m.pres_con_rev <- rma.mv(yi = partcorr, 
                         V = var_partcorr, 
                         mods = cbind(con.bstruc_rev, con.prop_rev, con.distm_rev, con.decpol_rev, con.decfis_rev),
                         slab = idstud,
                         data = edat.pres,
                         random = ~ 1 | idstud/idest, 
                         intercept = FALSE,
                         test = "t",
                         control=list(rel.tol=1e-8))
summary(m.pres_con_rev)
anova(m.pres_base, m.pres_con_rev, refit = TRUE)
fig5_1 <- m.pres_con %>% 
  tidy() %>% 
  bind_rows(tidy(m.pres_con_rev)) %>% 
  mutate(model = ifelse(grepl('_rev$', term), "Not included", "Included")) %>%
  mutate(term = str_remove(term, "_rev")) %>% 
  add_row(term = "Presidentialism", type = NA, estimate = NA, std.error = NA, statistic = NA, p.value = NA, model = "Included") %>% 
  add_row(term = "Presidentialism", type = NA, estimate = NA, std.error = NA, statistic = NA, p.value = NA, model = "Not included") %>% 
  relabel_predictors(c(con.bstruc = "Ballot Structure",
                       con.prop = "Electoral Formula",
                       con.distm = "District Magnitude",
                       con.decpol = "Pol. Decentralization",
                       con.decfis = "Fisc. Decentralization")) %>% 
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7),
         vars_order = c("Presidentialism", "Pol. Decentralization", "Fisc. Decentralization", "Electoral Formula", "District Magnitude", "Ballot Structure")) + 
  annotate("rect", xmin = 0.018, xmax = 0.138, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "Presidentialism") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  theme(legend.position = "none",
        plot.margin=grid::unit(c(0,1,0,0), "mm"),
        plot.title = element_text(face = "bold"),
        text = element_text(size = 15))

## COLUMN 2: POLITICAL DECENTRALIZATION
m.decpol_con <- rma.mv(yi = partcorr, 
                       V = var_partcorr, 
                       mods = cbind(con.bstruc, con.prop, con.distm, con.pres, con.decfis),
                       slab = idstud,
                       data = edat.decpol,
                       random = ~ 1 | idstud/idest,
                       intercept = FALSE,
                       test = "t",
                       control=list(rel.tol=1e-8))
summary(m.decpol_con)
anova(m.decpol_base, m.decpol_con, refit = TRUE)
m.decpol_con_rev <- rma.mv(y = partcorr, 
                           V = var_partcorr, 
                           mods = cbind(con.bstruc_rev, con.prop_rev, con.distm_rev, con.pres_rev, con.decfis_rev),
                           slab = idstud,
                           data = edat.decpol,
                           random = ~ 1 | idstud/idest, 
                           intercept = FALSE,
                           test = "t",
                           control=list(rel.tol=1e-8))
summary(m.decpol_con_rev)
anova(m.decpol_base, m.decpol_con_rev, refit = TRUE)
fig5_2 <- m.decpol_con %>% 
  tidy() %>% 
  bind_rows(tidy(m.decpol_con_rev)) %>% 
  mutate(model = ifelse(grepl('_rev$', term), "Not included", "Included")) %>%
  mutate(term = str_remove(term, "_rev")) %>% 
  add_row(term = "con.decpol", type = NA, estimate = NA, std.error = NA, statistic = NA, p.value = NA, model = "Included") %>% 
  add_row(term = "con.decpol", type = NA, estimate = NA, std.error = NA, statistic = NA, p.value = NA, model = "Not included") %>% 
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7),
         vars_order = c("con.pres", "con.decpol", "con.decfis", "con.prop", "con.distm", "con.bstruc")) + 
  annotate("rect", xmin = 0.013, xmax = 0.112, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "Pol. Decentralization") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  theme(axis.text.y=element_blank(),
        legend.position = "none",
        plot.margin=grid::unit(c(0,1,0,0), "mm"),
        plot.title = element_text(face = "bold"),
        text = element_text(size = 15))

## COLUMN 3: FISCAL DECENTRALIZATION
m.decfis_con <- rma.mv(y = partcorr, 
                       V = var_partcorr, 
                       mods = cbind(con.prop, con.pres, con.decpol),
                       slab = idstud,
                       data = edat.decfis,
                       random = ~ 1 | idstud/idest, 
                       intercept = FALSE,
                       test = "t",
                       control=list(rel.tol=1e-8))
summary(m.decfis_con)
anova(m.decfis_base, m.decfis_con, refit = TRUE)
m.decfis_con_rev <- rma.mv(y = partcorr, 
                           V = var_partcorr, 
                           mods = cbind(con.prop_rev, con.pres_rev, con.decpol_rev),
                           slab = idstud,
                           data = edat.decfis,
                           random = ~ 1 | idstud/idest, 
                           intercept = FALSE,
                           test = "t",
                           control=list(rel.tol=1e-8))
summary(m.decfis_con_rev)
anova(m.decfis_base, m.decfis_con_rev, refit = TRUE)
fig5_3 <- m.decfis_con %>% 
  tidy() %>% 
  bind_rows(tidy(m.decfis_con_rev)) %>% 
  mutate(model = ifelse(grepl('_rev$', term), "Not included", "Included")) %>%
  mutate(term = str_remove(term, "_rev")) %>% 
  add_row(term = "con.decfis", type = NA, estimate = NA, std.error = NA, statistic = NA, p.value = NA, model = "Included") %>% 
  add_row(term = "con.decfis", type = NA, estimate = NA, std.error = NA, statistic = NA, p.value = NA, model = "Not included") %>% 
  add_row(term = "con.distm", type = NA, estimate = NA, std.error = NA, statistic = NA, p.value = NA, model = "Included") %>% 
  add_row(term = "con.distm", type = NA, estimate = NA, std.error = NA, statistic = NA, p.value = NA, model = "Not included") %>% 
  add_row(term = "con.bstruc", type = NA, estimate = NA, std.error = NA, statistic = NA, p.value = NA, model = "Included") %>% 
  add_row(term = "con.bstruc", type = NA, estimate = NA, std.error = NA, statistic = NA, p.value = NA, model = "Not included") %>% 
  dwplot(dot_args = list(size = 2.5),
         whisker_args = list(size = 1),
         vline = geom_vline(xintercept = 0, color = "grey60", linetype = 2, size = 0.7),
         vars_order = c("con.pres", "con.decpol", "con.decfis", "con.prop", "con.distm", "con.bstruc")) + 
  annotate("rect", xmin = -0.175, xmax = -0.038, ymin = -Inf, ymax = Inf, alpha = .2) +
  labs(x = "", y = "", title = "Fisc. Decentralization") +
  scale_x_continuous(limits = c(-.475,.475),
                     breaks = c(-.4,-.2,0,.2,.4)) +
  scale_color_grey() +
  theme_bw() +
  theme(axis.text.y=element_blank(),
        legend.position = "none",
        plot.margin=grid::unit(c(0,1,0,0), "mm"),
        plot.title = element_text(face = "bold"),
        text = element_text(size = 15))

legend_con <- get_legend(fig5_1 + theme(legend.position="bottom", legend.title=element_blank()))
fig5 <- ggarrange(ggarrange(fig5_1, fig5_2, fig5_3, ncol = 3, widths = c(1.2,0.8,0.8)), 
                  legend_con,
                  nrow = 2, heights = c(1,.1))
fig5 <- ggdraw(fig5) +
  theme(plot.background = element_rect(fill="white", color = NA))

################################################################################

################################# APPENDICES ###################################

################ HETEROGENEITY ESTIMATES, APPENDIX B, TABLE 2 ###################


# Calculate overall model 
m_base <- rma.mv(y = partcorr, V = var_partcorr,
                 random = ~ 1 | idstud/idest, data = edat.all,
                 method = "ML")
summary(m_base)

# Extract level 2 and level 3 variances
level2_var <- m_base$sigma2[1]
level3_var <- m_base$sigma2[2]

print(level2_var)
print(level3_var)

# Get confidence intervals
ci_var <- confint(m_base)
print(ci_var)

# Get total variance
total_var <- level3_var + level2_var + mean(edat.all$var_partcorr)

# Calculate I^2 for each level
I2_level2 <- (level2_var / total_var) * 100
I2_level3 <- (level3_var / total_var) * 100


# Extract confidence intervals for variance components
level2_var_ci <- ci_var[[1]]$random
level2_var_ci_lb <- level2_var_ci[1, "ci.lb"]
level2_var_ci_ub <- level2_var_ci[1, "ci.ub"]

level3_var_ci <- ci_var[[2]]$random
level3_var_ci_lb <- level3_var_ci[1, "ci.lb"]
level3_var_ci_ub <- level3_var_ci[1, "ci.ub"]

# Compute confidence intervals for I^2
I2_level2_lower <- (level2_var_ci_lb / total_var) * 100
I2_level2_upper <- (level2_var_ci_ub / total_var) * 100

I2_level3_lower <- (level3_var_ci_lb / total_var) * 100
I2_level3_upper <- (level3_var_ci_ub / total_var) * 100

# Print the results
cat("I2_level2:", I2_level2, "\n95% CI:", I2_level2_lower, "-", I2_level2_upper, "\n")
cat("I2_level3:", I2_level3, "\n95% CI:", I2_level3_lower, "-", I2_level3_upper, "\n")

################################################################################

################ FIXED EFFECTS ANALYSIS, APPENDIX B, TABLE 3 ###################

## PRESIDENTIALISM
m.pres_base_fe <- rma.mv(yi = partcorr, 
                         V = var_partcorr,
                         slab = idstud,
                         data = edat.pres,
                         mods = ~ idstud, 
                         method = "FE", 
                         test = "t")
summary(m.pres_base_fe)
# Extract residual variance: 
residual_variance_pres <- vcov(m.pres_base_fe)[[1]]
# Calculate the total variance:
total_variance_pres <- sum(m.pres_base_fe$sigma2)
# Calculate the proportion of variability explained by the fixed effect (R^2)
R2_value_pres <- 1 - (residual_variance_pres / total_variance_pres)

## POLITICAL DECENTRALIZATION
m.decpol_base_fe <- rma.mv(yi = partcorr, 
                           V = var_partcorr,
                           slab = idstud,
                           data = edat.decpol,
                           mods = ~ idstud, 
                           method = "FE", 
                           test = "t")
summary(m.decpol_base_fe)
residual_variance_decpol <- vcov(m.decpol_base_fe)[[1]]
total_variance_decpol <- sum(m.decpol_base_fe$sigma2)
R2_value_decpol <- 1 - (residual_variance_decpol / total_variance_decpol)

## FISCAL DECENTRALIZATION
m.decfis_base_fe <- rma.mv(yi = partcorr, 
                           V = var_partcorr,
                           slab = idstud,
                           data = edat.decfis,
                           mods = ~ idstud, 
                           method = "FE", 
                           test = "t")
summary(m.decfis_base_fe)
residual_variance_decfis <- vcov(m.decfis_base_fe)[[1]]
total_variance_decfis <- sum(m.decfis_base_fe$sigma2)
R2_value_decfis <- 1 - (residual_variance_decfis / total_variance_decfis)

################################################################################

########### WITHIN AND ACROSS STUDY VARIATION, APPENDIX B, TABLE 4 #############

## PRESIDENTIALISM
tau_sq <- vcov(m.pres_base)[[1]]
# Extract the total observed variance
total_variance <- sum(m.pres_base$sigma2)
# Calculate the proportion of total variance due to within-study variability (I^2_2)
I2_within <- (tau_sq / total_variance) * 100
# Calculate the proportion of total variance due to between-study variability (I^2_3)
I2_between <- 100 - I2_within
# Print the results
cat("Within-Study Variation (I^2_2):", I2_within, "%\n")
cat("Between-Study Variation (I^2_3):", I2_between, "%\n")

## POLITICAL DECENTRALIZATION
tau_sq <- vcov(m.decpol_base)[[1]]
total_variance <- sum(m.decpol_base$sigma2)
I2_within <- (tau_sq / total_variance) * 100
I2_between <- 100 - I2_within
cat("Within-Study Variation (I^2_2):", I2_within, "%\n")
cat("Between-Study Variation (I^2_3):", I2_between, "%\n")

## FISCAL DECENTRALIZATION
tau_sq <- vcov(m.decfis_base)[[1]]
total_variance <- sum(m.decfis_base$sigma2)
I2_within <- (tau_sq / total_variance) * 100
I2_between <- 100 - I2_within
cat("Within-Study Variation (I^2_2):", I2_within, "%\n")
cat("Between-Study Variation (I^2_3):", I2_between, "%\n")

################################################################################

####### NUMERICAL RESULTS FOR HETEROGENOUS EFFECTS APPENDIX C, TABLE 7 #########

## PRESIDENTIALISM
## PUBLISHED
m.pres_pub <- rma.mv(yi = partcorr, 
                     V = var_partcorr, 
                     mods = cbind(pub.yes, pub.no),
                     slab = idstud,
                     data = edat.pres,
                     random = ~ 1 | idstud/idest, 
                     intercept = FALSE,
                     test = "t")
summary(m.pres_pub)
anova(m.pres_base, m.pres_pub, refit = TRUE)
## TOP 10 JOURNAL
m.pres_top10 <- rma.mv(yi = partcorr, 
                       V = var_partcorr, 
                       mods = cbind(top10.yes, top10.no),
                       slab = idstud,
                       data = edat.pres,
                       random = ~ 1 | idstud/idest, 
                       intercept = FALSE,
                       test = "t")
summary(m.pres_top10)
anova(m.pres_base, m.pres_top10, refit = TRUE)
## DISCIPLINE
m.pres_field <- rma.mv(yi = partcorr, 
                       V = var_partcorr, 
                       mods = cbind(field.econ, field.polsci, field.other, field.multi),
                       slab = idstud,
                       data = edat.pres,
                       random = ~ 1 | idstud/idest, 
                       intercept = FALSE,
                       test = "t")
summary(m.pres_field)
anova(m.pres_base, m.pres_field, refit = TRUE)
## PUBLICATION YEAR
m.pres_ypub <- rma.mv(yi = partcorr, 
                      V = var_partcorr, 
                      mods = cbind(ypub.99_04, ypub.05_09, ypub.10_14, ypub.15_19, ypub.20_24),
                      slab = idstud,
                      data = edat.pres,
                      random = ~ 1 | idstud/idest, 
                      intercept = FALSE,
                      test = "t")
summary(m.pres_ypub)
anova(m.pres_base, m.pres_ypub, refit = TRUE)
## NO. AUTHORS
m.pres_autno <- rma.mv(yi = partcorr, 
                       V = var_partcorr, 
                       mods = cbind(auth_no.1, auth_no.2, auth_no.3, auth_no.4),
                       slab = idstud,
                       data = edat.pres,
                       random = ~ 1 | idstud/idest, 
                       intercept = FALSE,
                       test = "t")
summary(m.pres_autno)
anova(m.pres_base, m.pres_autno, refit = TRUE)

####### NUMERICAL RESULTS FOR HETEROGENOUS EFFECTS APPENDIX C, TABLE 8 #########

## POLITICAL DECENTRALIZATION
m.decpol_base_comp <- rma.mv(yi = partcorr, 
                             V = var_partcorr,
                             slab = idstud,
                             data = edat.decpol %>% filter(!is.na(pub.yes)),
                             random = ~ 1 | idstud/idest, 
                             method="ML")
summary(m.decpol_base_comp)
## PUBLISHED
m.decpol_pub <- rma.mv(yi = partcorr, 
                       V = var_partcorr, 
                       mods = cbind(pub.yes, pub.no),
                       slab = idstud,
                       data = edat.decpol %>% filter(!is.na(pub.yes)),
                       random = ~ 1 | idstud/idest, 
                       intercept = FALSE,
                       method = "ML")
summary(m.decpol_pub)
anova(m.decpol_base_comp, m.decpol_pub, refit = TRUE)
## TOP 10 JOURNAL
m.decpol_top10 <- rma.mv(yi = partcorr, 
                         V = var_partcorr, 
                         mods = cbind(top10.yes, top10.no),
                         slab = idstud,
                         data = edat.decpol %>% filter(!is.na(pub.yes)),
                         random = ~ 1 | idstud/idest, 
                         intercept = FALSE,
                         method = "ML")
summary(m.decpol_top10)
anova(m.decpol_base_comp, m.decpol_top10, refit = TRUE)
## DISCIPLINE
m.decpol_field <- rma.mv(yi = partcorr, 
                         V = var_partcorr, 
                         mods = cbind(field.econ, field.polsci, field.other, field.multi),
                         slab = idstud,
                         data = edat.decpol %>% filter(!is.na(pub.yes)),
                         random = ~ 1 | idstud/idest, 
                         intercept = FALSE,
                         method = "ML")
summary(m.decpol_field)
anova(m.decpol_base_comp, m.decpol_field, refit = TRUE)
## PUBLICATION YEAR
m.decpol_ypub <- rma.mv(yi = partcorr, 
                        V = var_partcorr, 
                        mods = cbind(ypub.99_04, ypub.05_09, ypub.10_14, ypub.15_19, ypub.20_24),
                        slab = idstud,
                        data = edat.decpol %>% filter(!is.na(pub.yes)),
                        random = ~ 1 | idstud/idest, 
                        intercept = FALSE,
                        method = "ML")
summary(m.decpol_ypub)
anova(m.decpol_base_comp, m.decpol_ypub, refit = TRUE)
## NO. AUTHORS
m.decpol_autno <- rma.mv(yi = partcorr, 
                         V = var_partcorr, 
                         mods = cbind(auth_no.1, auth_no.2, auth_no.3, auth_no.4),
                         slab = idstud,
                         data = edat.decpol,
                         random = ~ 1 | idstud/idest, 
                         intercept = FALSE,
                         method = "ML")
summary(m.decpol_autno)
anova(m.decpol_base_comp, m.decpol_autno, refit = TRUE)

####### NUMERICAL RESULTS FOR HETEROGENOUS EFFECTS APPENDIX C, TABLE 9 #########

## FISCAL DECENTRALISATION
m.decfis_base_comp <- rma.mv(yi = partcorr, 
                             V = var_partcorr,
                             slab = idstud,
                             data = edat.decfis %>% filter(!is.na(pub.yes)),
                             random = ~ 1 | idstud/idest, 
                             method="ML")
summary(m.decfis_base_comp)
## PUBLISHED
m.decfis_pub <- rma.mv(yi = partcorr, 
                       V = var_partcorr, 
                       mods = cbind(pub.yes, pub.no),
                       slab = idstud,
                       data = edat.decfis %>% filter(!is.na(pub.yes)),
                       random = ~ 1 | idstud/idest, 
                       intercept = FALSE,
                       method = "ML")
summary(m.decfis_pub)
anova(m.decfis_base_comp, m.decfis_pub, refit = TRUE)
## TOP 10 JOURNAL
m.decfis_top10 <- rma.mv(yi = partcorr, 
                         V = var_partcorr, 
                         mods = cbind(top10.yes, top10.no),
                         slab = idstud,
                         data = edat.decfis %>% filter(!is.na(pub.yes)),
                         random = ~ 1 | idstud/idest, 
                         intercept = FALSE,
                         method = "ML")
summary(m.decfis_top10)
anova(m.decfis_base_comp, m.decfis_top10, refit = TRUE)
## DISCIPLINE
m.decfis_field <- rma.mv(yi = partcorr, 
                         V = var_partcorr, 
                         mods = cbind(field.econ, field.polsci, field.other, field.multi),
                         slab = idstud,
                         data = edat.decfis %>% filter(!is.na(pub.yes)),
                         random = ~ 1 | idstud/idest, 
                         intercept = FALSE,
                         method = "ML")
summary(m.decfis_field)
anova(m.decfis_base_comp, m.decfis_field, refit = TRUE)
## PUBLICATION YEAR
m.decfis_ypub <- rma.mv(yi = partcorr, 
                        V = var_partcorr, 
                        mods = cbind(ypub.99_04, ypub.05_09, ypub.10_14, ypub.15_19, ypub.20_24),
                        slab = idstud,
                        data = edat.decfis %>% filter(!is.na(pub.yes)),
                        random = ~ 1 | idstud/idest, 
                        intercept = FALSE,
                        method = "ML")
summary(m.decfis_ypub)
anova(m.decfis_base_comp, m.decfis_ypub, refit = TRUE)
## NO. AUTHORS
m.decfis_autno <- rma.mv(yi = partcorr, 
                         V = var_partcorr, 
                         mods = cbind(auth_no.1, auth_no.2, auth_no.3, auth_no.4),
                         slab = idstud,
                         data = edat.decfis %>% filter(!is.na(pub.yes)),
                         random = ~ 1 | idstud/idest, 
                         intercept = FALSE,
                         method = "ML")
summary(m.decfis_autno)
anova(m.decfis_base_comp, m.decfis_autno, refit = TRUE)

################################################################################